function A305()
    format long;
    % 谱估计根据随机过程的有限观测值确定该过程的谱，包括：经典谱估计、参数建模、最大熵谱估计等。

    % 读取Excel文件
    data_single=readtable('A3_data_single.xlsx');
    data_double =readtable('A3_data_double.xlsx');
    xn1 =data_single{:,150};              %读取第150列(单排钢筋衬砌信号)
    xn2 =data_double{:,150};              %读取第150列(双排钢筋衬砌信号)
    
    % 谱估计参数
    Fs = 1000; 
     % 图1：单排钢筋信号处理与分析
    figure(1);
    process_signal(xn1, xn2, Fs, '单排钢筋信号', '双排钢筋信号');
end
    
function process_signal(xn1, xn2, Fs,signal_name1, signal_name2)
    % 信号时域生成
    N1 = length(xn1);  % 获取信号1的总长度
    N2 = length(xn2);  % 获取信号2的总长度
    Nsec = 128; Nfft = Nsec; 
    
    % 确保信号分段一致，取两信号的最小长度
    N = min(N1, N2);  % 确保信号1和信号2一致处理
    num_segments = floor(N / Nsec);  % 计算可以分为多少段
    
    % 时域信号显示
    n1 = 0:N1-1; t1 = n1/Fs;
    n2 = 0:N2-1; t2 = n2/Fs;
    subplot(3, 3, 1);
    plot(t1, xn1, 'k-', 'DisplayName', signal_name1); hold on;
    plot(t2, xn2, 'r-', 'DisplayName', signal_name2); hold off;
    xlabel('时间 (秒)'); ylabel('信号值'); 
    title('时域信号');
    legend;

    % 周期图法
    P1 = 10*log10(abs(fft(xn1(1:Nfft), Nfft).^2)/Nfft);
    P2 = 10*log10(abs(fft(xn2(1:Nfft), Nfft).^2)/Nfft);
    f = (0:length(P1)-1) * Fs/length(P1);
    subplot(3, 3, 2);
    plot(f(1:Nsec/2), P1(1:Nsec/2), 'k-', 'DisplayName', signal_name1); hold on;
    plot(f(1:Nsec/2), P2(1:Nsec/2), 'r-', 'DisplayName', signal_name2); hold off;
    title('周期图 N=128（dB图）');
    legend;
    ylabel('功率谱/dB');
    xlabel('频率/Hz');

   % 平均周期图法
    Pxx1 = zeros(1, Nsec);  % 初始化信号1的功率谱容器
    Pxx2 = zeros(1, Nsec);  % 初始化信号2的功率谱容器

% 遍历信号1每段数据进行FFT计算
    for i = 1:num_segments
        segment_start = (i - 1) * Nsec + 1;  % 当前段的起始位置
        segment_end = i * Nsec;  % 当前段的结束位置
        segment_data1 = xn1(segment_start:segment_end);  % 提取信号1的当前段数据
        segment_data2 = xn2(segment_start:segment_end);  % 提取信号2的当前段数据
    
        % 如果最后一段不足128个点，可以选择零填充
        if length(segment_data1) < Nsec
            segment_data1 = [segment_data1, zeros(1, Nsec - length(segment_data1))];
        end
        if length(segment_data2) < Nsec
            segment_data2 = [segment_data2, zeros(1, Nsec - length(segment_data2))];
        end
    
        % 计算当前段的功率谱
        Pxx_segment1 = abs(fft(segment_data1, Nsec).^2) / Nsec;  % 计算信号1每段的功率谱
        Pxx_segment2 = abs(fft(segment_data2, Nsec).^2) / Nsec;  % 计算信号2每段的功率谱
    
        % 累加功率谱
        Pxx1 = Pxx1 + Pxx_segment1;
        Pxx2 = Pxx2 + Pxx_segment2;
    end

    % 求取平均功率谱并转化为dB
    Pxx1 = 10 * log10(Pxx1 / num_segments);  % 转化为dB并取平均
    Pxx2 = 10 * log10(Pxx2 / num_segments);  % 转化为dB并取平均

    % 计算频率轴
    f = (0:length(Pxx1)-1) * Fs / length(Pxx1);  % 给出频率序列

    % 绘制平均功率谱（两个信号在同一图）
    subplot(3, 3, 3);
    plot(f(1:Nsec/2), Pxx1(1:Nsec/2), 'k-', 'DisplayName', signal_name1); hold on;
    plot(f(1:Nsec/2), Pxx2(1:Nsec/2), 'r-', 'DisplayName', signal_name2); hold off;
    legend;
    xlabel('频率/Hz');
    ylabel('功率谱/dB');
    title(['平均周期图法 Nsec=' num2str(Nsec)]);

% 修改后的平均周期图法
% Hamming窗口
M = 128;  % 每段长度
window = hamming(M);  % 128点的Hamming窗口
Pxxtotal1 = 0;
Pxxtotal2 = 0;
L = floor(N / (M / 2)) - 1;  % 计算段数，采用重叠段方式
for i = 1:L-1
    % 当前段的索引（根据重叠段调整）
    segment_start = (i - 1) * (M / 2) + 1;  % 每段起始位置
    segment_end = segment_start + M - 1;  % 每段结束位置
    segment_data1 = xn1(segment_start:segment_end);  % 提取当前段数据
    segment_data2 = xn2(segment_start:segment_end);  % 提取当前段数据
    m1 = abs(fft(window' .* segment_data1, M).^2) / M;
    m2 = abs(fft(window' .* segment_data2, M).^2) / M;
    Pxxtotal1 = Pxxtotal1 + m1;  % 累加每段的功率谱
    Pxxtotal2 = Pxxtotal2 + m2;  % 累加每段的功率谱
end
% 最后一段处理
window_last = hamming(M);  % 最后一段的窗口
segment_last_start = N - M + 1;  % 最后一段的起始位置
mend1 = abs(fft(window_last' .* xn1(segment_last_start:N), M).^2) / M;  % 处理最后一段
mend2 = abs(fft(window_last' .* xn2(segment_last_start:N), M).^2) / M;  % 处理最后一段
Pxxtotal1 = (Pxxtotal1 + mend1) / L;  % 计算所有段的平均功率谱
Pxxtotal2 = (Pxxtotal2 + mend2) / L;  % 计算所有段的平均功率谱
% 转化为dB
Pxx1 = 10*log10(Pxxtotal1);
Pxx2 = 10*log10(Pxxtotal2);
% 计算频率序列
f1 = (0:length(Pxx1)-1) * Fs / length(Pxx1);
f2 = (0:length(Pxx2)-1) * Fs / length(Pxx2);
% 给出频率序列
w1 = fft(window,N);
w10 = abs(fftshift(w1));
% 绘制功率谱曲线
subplot(3, 3, 4);
plot(f1(1:M/2), Pxx1(1:M/2), 'k-', 'DisplayName', signal_name1);hold on;  % 绘制功率谱曲线
plot(f2(1:M/2), Pxx2(1:M/2), 'r-', 'DisplayName', signal_name1);hold off;  % 绘制功率谱曲线
title('Hamming');
legend;
xlabel('频率/Hz'); ylabel('功率谱/dB');

    
    

% 相关功率谱估计（BT法）
nfft = 128;
ncxk = 3*nfft/4;
cxn1 = xcorr(xn1, 'unbiased');
cxn2 = xcorr(xn2, 'unbiased');
CXk1 = fft(cxn1, ncxk);
CXk2 = fft(cxn2, ncxk);
Pxx1 = 10*log10(abs(CXk1));
Pxx2 = 10*log10(abs(CXk2));
index = 0:round(ncxk/2-1);
k = index*Fs/ncxk;
C1 = Pxx1(index+1);
C2 = Pxx1(index+1);
P1 = (Pxx1(index+1));
P2 = (Pxx2(index+1));
subplot(3, 3, 5)
plot(k, P1, 'k-', 'DisplayName', signal_name1);hold on;
plot(k, P2, 'r-', 'DisplayName', signal_name2);hold off;
%var(C)
title('BT法功率谱估计 N=128');
legend;
xlabel('频率 Hz'); ylabel('功率谱/dB');

% AR模型参数估计
ys = xn1;  
yd = xn2; 
ybars = mean(ys); 
ybard = mean(yd); 
zs = ys - ybars;  % 将信号去均值
zd = yd - ybard;  % 将信号去均值
% 对1-15阶AR模型逐个参数估计，在工作区m/Report/Fit，分别进行以下参数估计
y1_s = []; y2_s = []; y3_s = [];  % 信号1的FPE、AIC、BIC
y1_d = []; y2_d = []; y3_d = [];  % 信号2的FPE、AIC、BIC
% 采用循环分别为y1、y2和y3赋值
for i = 1:1:15
    % 对信号1进行AR模型估计
    m1 = ar(zs, i, 'fb');  
    y1_s = [y1_s m1.Report.Fit.FPE];
    y2_s = [y2_s m1.Report.Fit.AIC];
    y3_s = [y3_s m1.Report.Fit.BIC];
    
    % 对信号2进行AR模型估计
    m2 = ar(zd, i, 'fb');
    y1_d = [y1_d m2.Report.Fit.FPE];
    y2_d = [y2_d m2.Report.Fit.AIC];
    y3_d = [y3_d m2.Report.Fit.BIC];
end

% 绘制 AR 模型阶数选择图
x = (1:1:15);
% 清空子图，避免累积
subplot(3, 3, 6);
cla;
% 绘制信号1和信号2的结果
subplot(3, 3, 6);
yyaxis left
plot(x, y1_s, 'b*-', 'DisplayName', '单排钢筋-FPE'); hold on;
plot(x, y1_d, 'g*-', 'DisplayName', '双排钢筋-FPE');
ylabel('FPE');  % 左轴为FPE
yyaxis right
plot(x, y2_s, 'ro-', 'DisplayName', '单排钢筋-AIC');
plot(x, y3_s, 'rs-', 'DisplayName', '单排钢筋-BIC');
plot(x, y2_d, 'mo-', 'DisplayName', '双排钢筋-AIC');
plot(x, y3_d, 'ys-', 'DisplayName', '双排钢筋-BIC');hold off;
ylabel('AIC/BIC');  % 右轴为AIC/BIC
xlabel('AR模型阶数');
legend;
title('单排钢筋和双排钢筋的AR模型参数估计');

% AR(p)过程的Burg PSD估计
% 使用 Burg 方法进行 AR(p) 模型参数估计
p = 6; % 模型阶数
% 对信号1和信号2分别估计AR(p)模型
[A1, ~] = arburg(zs, p);  % 估计信号1的AR(p)模型参数
[A2, ~] = arburg(zd, p);  % 估计信号2的AR(p)模型参数

% 使用pburg估计信号的功率谱
[Pxx1, F_pburg1] = pburg(xn1, 6, length(xn1), 1);  % 使用pburg估计信号1的功率谱
[Pxx2, F_pburg2] = pburg(xn2, 6, length(xn2), 1);  % 使用pburg估计信号2的功率谱
% 绘制AR过程的频率响应（freqz）
subplot(3, 3, 7);
[H1, F1] = freqz(1, A1, [], 1);  % 计算信号1的频率响应
[H2, F2] = freqz(1, A2, [], 1);  % 计算信号2的频率响应
% 绘制频率响应
plot(F1, 10*log10(abs(H1)), 'b-', 'DisplayName', '单排钢筋的功率谱密度'); hold on;
plot(F2, 10*log10(abs(H2)), 'r-', 'DisplayName', '双排钢筋的功率谱密度');
% 绘制功率谱估计
plot(F_pburg1, 10*log10(Pxx1), 'k-', 'DisplayName', '单排钢筋的pburg PSD估计'); 
plot(F_pburg2, 10*log10(Pxx2), 'g-', 'DisplayName', '双排钢筋的pburg PSD估计');
% 添加图例和标签
legend;
xlabel('频率（Hz）');
ylabel('功率谱 PSD（dB/Hz）');
title('AR过程的Burg PSD估计');
hold off;

% 用于模型阶次确定的反射系数
rng default  % 设置相同的随机数种子
x1 = filter(1, A1, randn(1000, 1));  % 使用 AR(p) 模型生成随机信号1
x2 = filter(1, A2, randn(1000, 1));  % 使用 AR(p) 模型生成随机信号2
% 对信号1进行AR(12)模型反射系数估计
[a1, e1, k1] = arburg(x1, 12);  
% 对信号2进行AR(12)模型反射系数估计
[a2, e2, k2] = arburg(x2, 12);
subplot(3, 3, 8)
stem(k1, 'filled')  % 绘制信号1的反射系数
hold on;
stem(k2, 'filled')  % 绘制信号2的反射系数
hold off;
title('反射系数')
xlabel('模型阶次');
ylabel('采样自相关');
legend({'信号1的反射系数', '信号2的反射系数'});

% 多通道信号的Burg PSD估计
Fs = 1000;  % 采样频率
t = 0:1/Fs:1-1/Fs;  % 生成1秒的时间向量
% 使用AR(p)模型生成信号
x1 = filter(1, A1, randn(length(t), 1));  % 使用 AR(p) 模型生成信号1
x2 = filter(1, A2, randn(length(t), 1));  % 使用 AR(p) 模型生成信号2
% 生成一个新的多通道信号（这里将x1和x2作为两个通道信号）
% 加上一些噪声和不同频率的成分
f = [100; 200; 300];  % 频率成分
x3 = cos(2*pi*f(1)*t)';  % 生成第三个信号
x_multi = [x1, x2, x3] + randn(length(t), 3);  % 合并信号形成多通道信号
% 使用pburg对多通道信号进行功率谱估计
morder = 12;  % 模型阶数
subplot(3, 3, 9)
pburg(x_multi, morder, [], Fs)
legend('单排钢筋信号', '双排钢筋信号', '合成信号 ')
title('多通道信号的Burg PSD估计')
end
